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o : 

. We explore, in three spatial dimensions, the transition from the normal state to the Fulde-Ferrel- 

£\j ■ Larkin-Ovchinnikov superfluid phases. We restrict ourselves to the case of the 'planar' phase, where 

the order parameter depends only on a single spatial coordinate. We first show that, in the case 
of the simple Fulde-Ferrell phase, singularities occur at zero temperature in the free energy which 
prevents, at low temperature, a reliable use of an expansion in powers of the order parameter. We 
then introduce in the quasiclassical equations a Fourier expansion for the order parameter and the 
Green's functions, and we show that it converges quite rapidly to the exact solution. We finally 
implement numerically this method and find results in excellent agreement with the earlier work of 
Matsuo et al. In particular when the temperature is lowered from the tricritical point, the transition 
switches from first to second order. In the case of the first order transition, the spatial dependence 
of the order parameter at the transition is found to be always very nearly a pure cosine, although 
the maximum of its modulus may be comparable to the one of the uniform BCS phase. 

CLf PACS numbers : 74.20.Fg, 74.60. Ec 
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I. INTRODUCTION 



Despite being actively investigated for fourty years the problem of the structure of the superconducting order 
Ch ' parameter in very high fields is still the subject of intensive research. In the compounds of interest the coupling of 
the magnetic field to electronic spins can no longer be ignored, and the situation where it is the only relevant one 
i^i - has to be considered. In this case one faces the problem of pairing electrons, for which the spin up and spin down 
chemical potentials are not the same. This question has been addressed independently by Fulde and Ferrell [1] (FF) 
t-H , and by Larkin and Ovchinnikov [2] (LO), who proposed that the best order parameter corresponds to pairs formed 
with a non zero total momentum, in contrast to the standard situation of the BCS theory. It is worth noting that 
this kind of problem has been found recently to be quite relevant for ultracold atomic Fermi gases [3] as well as for 
the physics of neutron stars [4-6]. More specifically Larkin and Ovchinnikov [2] considered for the order parameter 
superpositions of different plane waves, corresponding physically to different pair total momentum. They investigated 
_ ■ which superposition was favored near the transition at T = 0. Nevertheless they considered only a second order phase 
transition, which is not the most general situation as we will discuss below. Hence the question of the exact structure 
of the order parameter in the FFLO phases is still an open problem. Since most experiments identifying tentatively 
FFLO phases rest heavily on the theoretical analysis, this is also a problem with major experimental implications. 

In a preceding paper [7] we have investigated analytically the transition to the FFLO phases in the vicinity of the 
tricritical point (TCP), where the FFLO transition line starts. This point is located at T tcp /T c o = 0.561 where T c o 
is the critical temperature for fi, = , with 2/2 = ^ — /ij being the chemical potential difference between the two 
fermionic populations forming pairs, as for example spin up and down electrons (the corresponding effective field is 
prtcp/T c o — 1.073). In agreement with preceding numerical work [8-10] we have found that the transition is first order 
to an order parameter which is, to a very good precision, simply proportional to a one-dimensional 'planar' texture 
A(r) ~ cos(q.r). This order parameter is actually the one which has been shown by Larkin and Ovchinnikov [2] to 
• i— i . be the most favorable for a (second order) transition at T = 0. Compared with other works we have been able to 
understand qualitatively and quantitatively the reasons which favor this order parameter with respect to all the other 
possible ones. Namely we have shown that a real order parameter is favored, and that, among these states, those with 
the smallest number of plane waves are prefered. This then leads to an order parameter with a cos(qo.r) dependence, 
in agreement with preceding work. 

A remarkable feature of the results is that the location in the p,, T plane of this first order transition toward the 
'planar' order parameter is very near the standard FFLO second order phase transition. This is true not only near the 
TCP [7-9] but also almost down to T = [10] (actually the transition to the 'planar' order parameter goes back to a 
second order phase transition at low temperature in agreement with LO). This proximity of a second order transition 
may lead to believe that the order parameter A is reasonably small at the first order transition. This is trivially valid 
near the TCP where a Landau- Ginzburg type expansion up to sixth order in order parameter A could be performed 
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[7,8], but it is a tempting hypothesis even at lower temperature. This possibility has been explored by Houzet et 
al. [9]. They found problems in applying this scheme because, for the stablest phase, namely the planar one, the 
coefficient of the sixth order in A changes sign when the temperature is lowered not much below the TCP. This leads 
to an instability and so the expansion up to sixth order in powers of A becomes inconsistent. 

We will first analyze this problem and show that it is already present in the simple case of the Fulde-Ferrell phase 
where the order parameter is given by A(r) = A exp(iq.r). This analysis gives a clear hint that an expansion in powers 
of A is going to fail anyway at low temperature. This suggest that one should avoid to perform such an expansion. 
Among the various possibilities for improving the situation, one of them is to remark that, at the transition, near 
the TCP, the actual order parameter is quite near a simple superposition of plane waves [7] even if the transition is 
first order (for a second order transition the order parameter is exactly such a superposition, as investigated [2] for 
example by Larkin and Ovchinnikov at T = 0). So the power expansion near the TCP amounts also to keep only the 
lowest order in a Fourier expansion of the order parameter. This leads to look for a Fourier expansion in the equations 
instead of a A expansion. 

Since we want to deal with the full nonlinear, space-dependent problem, the convenient starting point is not 
Gorkov's equations, but rather the quasiclassical equations of Eilenberger [11], Larkin and Ovchinnikov [12]. Not only 
are those equations in their simplest form the most compact and convenient formulation of our problem, but a major 
advantage is that they can be extended in full generality to much more complex situations [13] and allow to formulate 
transport problems, including many-body effects, with the same level of efficiency However since we have to deal 
with a comparatively simpler problem, we will use for simplicity in this paper the original formulation and notations 
of Eilenberger [11]. In comparison the general formalism is used by Burkardt and Raincr [14] for an analysis of a 
FFLO transition in two dimensions with a planar order parameter. 

In this paper we will show that the introduction of a Fourier expansion in the quasiclassical equations allows to 
obtain a solution which converges very rapidly toward the exact result. As a consequence a few terms in the expansion 
provide an excellent approximation. Here we will just deal with the principle of this method and its application to the 
planar transition. In particular we will rederive the results of Matsuo et al [10]. Applications to other more complex 
cases, which are the more fundamental interest of this procedure, will be considered in another paper. 

The paper is organized as follows. In the next section we consider the free energy and study in particular the simple 
case of the Fulde-Ferrell phase, and show that it has singularities at T — which make an expansion in powers of 
the order parameter unreliable. In section III we explain our Fourier expansion for the simplest case of a cosine order 
parameter. This is then generalized in the following section to any (one-dimensional) order parameter. Finally we 
give in section V the results of the numerical implementation of our method. 



II. THE FULDE-FERRELL PHASE T = FREE ENERGY 



We will show that the problems arising in the expansion of the free energy in powers of the order parameter are 
already present when one considers the simple Fulde-Ferrell (FF) state. Let us start with the completely general 
expression, that we will use further on, for the free energy difference per unit volume between the superconducting 
state and the normal state [14-16] : 

dr-\A(r)\ 2 + 4ttTN RcJ2 du / — ±\g a (u>, k, r) - g n (u, k, r)] (1) 

Here V is the standard BCS interaction and Nq is the single spin density of states at the Fermi energy. The difference 
/if — Hi = 2/2 between spin up and spin down chemical potentials comes in the definition of u> n = w n — ifl where 
uj n = ttT(2ti + 1) are Matsubara frequencies. For the " £- integrated " or quasiclassical Green's functions we have 
used Eilcnberger's notations g(u,)s., r) = — J d£kG(u>, k, r) where £k is the kinetic energy measured from the average 
Fermi level (1/2) (/if + and G(ui n , k, r) is the usual temperature Green's function (these Green's functions we deal 
with are those for up spin, the down spin Green's functions are obtained by a simple transform and the sum over 
the spin leads to take the real part in Eq.(l)). With these notations we have g n (u> n ,i.,r) = 1 for cu n > 0. It results 
directly from the starting Gorkovs equations [2] that the Green's functions in the presence of the effective field ft are 
obtained from those in the absence of it by the simple replacement of w n by w n . 

We look now for the expression of this free energy at T = for the simple FF state where the order parameter is 
given by A(r) = Aexp(iq.r). Since in this case |A(r)| 2 = A 2 , the Green's function is just obtained from the standard 
BCS one by shifting [2] all the momenta by q/2. Finally the quasiclassical Green's function is just the BCS one, 
except that we have to change ui into ui — ijik with fik = m(1 — ?k.q), where we have defined the reduced wavevector 
q = gfcir/(2m/x) (this results also from the general Eilenberger's equations we will write below). 
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The free energy for the standard uniform BCS phase with \i = is : 

n = n s - n n = ^n a 2 x 2 \n(x 2 /e) (2) 

where Ao = 2uj c exp(— 1/NqV) is the zero temperature BCS phase gap (u c is the standard cut-off of BCS theory), 
and we have expressed A in units of A by introducing x = A/A . This free energy is naturally minimum for x = 1 
and the minimum is — ^N A 2 . In the presence of a nonzero effective field p > this expression becomes from Eq.(l): 

Q X 2 

rh 2 + Re[x 2 ln(m + V m 2 — x 2 ) — my m 2 — x 2 } (3) 



N A 2 2 

where we have also expressed p in units of A by setting rh = P/Aq. For A > p this free energy reduces to 
7 = ix 2 ln(x 2 /e) + rh 2 and gives the standard Clogston-Chandrasekhar [17,18] first order transition m = l/y/2. 
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On the other hand it gives for small A -C p the expansion M n A i = x 2 ln(2m) — x A l&fh 2 — x 6 /32m 4 , leading in particular 
to the second order spinodal transition for m = 1/2. This expansion can be generalized at T ^ as : 

§- = ln[T/T sp (p/T)} A 2 + E(-1F + A 2p A 2 ?+ 2 (4) 

with : 

i 

A 2p = 2nTBe[J2^2 P +i} (5) 

n=0 W « 

and T sp (p/T) is the temperature of the second order spinodal transition toward the standard BCS phase. It is 
interesting to note that, while the coefficients A 2p are clearly all positive when p — > 0, they are given by A 2p = 
(— l) p /(2pp 2p ) when T — > 0. Moreover one can see that A 2p has p zeros when p/T goes from to oo (one goes 
basically from A 2p to A 2p+2 by taking a double derivative with respect to p) . Hence the higher order coefficients have 
many changes of sign in the low temperature range. This feature corresponds to the singular behaviour which occurs 
for A = p at T = in Eq.(3). It allows also to understand that the change of signs found by Houzet et al. [9] arc not 
simple accidents, but a systematic behaviour linked to the singularity appearing at T = 0. 

Finally the T = free energy of the FF phase is obtained by replacing in Eq.(3) p by pk = p{l — ?k.q) and 
averaging over the direction k as in Eq.(l). We give the result only in the case where q > 1 since this is the range of 
wavevector corresponding to the actual minimum of the free energy. One finds : 
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— + m 2 (l + — ) + — ^Rc[m + ln(m + + y fa\ — x 2 ) — \jfh\ — x 2 + (rh + 
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•^Re[(m 2 -x 2 f' 2 + {fh + -»m_)l (6) 
6mq 



where we have used the notation rh± = rh(q ± 1). 

This result has singularities for A = p(q ± 1). These are just the manifestation of the singularity found in Eq.(3) 
for A = p, corresponding to the upper and lower bound in the angular integration. In particular the singularity at 
A = p(q — 1) gives the radius of convergence of the expansion in powers of A. A particular consequence is that no 
expansion is possible for q = 1. This is just the situation which is found when one works in a two-dimensional space. 
This singular situation leads to a cascade involving an infinite number of phase transitions when the temperature goes 
to zero, as we have shown elsewhere [19]. In the case of a three-dimensional space, with which we deal in this paper, 
the radius of convergence is non zero, but it is fairly small since the minimum free energy is found at low temperature 
for values of q not far from the T = LO result q = 1.2. Therefore a rapidly convergent expansion for the free energy 
is only valid for quite small A, and this happens to be in contradiction with the values of A needed to minimize this 
free energy. Naturally this expansion of Eq. (6) can be performed explicitely and the problem with the convergence is 
then quite obvious. 

Now it is clear that these same problems arise if, instead of a phase with a single plane-wave as is the FF phase, 
we consider a more complicated phase which is a sum of plane waves, such as the planar phase A(r) ~ cos(q.r). This 
is already obvious from the fact that the terms which arise in the expansion for the FF phase will also appear in the 
expansion for this phase. Other terms with weaker singularities at A = p(q — 1) will also be present. We note that 
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a singularity is already present in the fourth order terms investigated [2] by LO, as it can be seen from the explicit 
expression of their integral J, but it occurs for a specific value of the angle between the wavevectors which happens to 
be irrelevant for their final conclusion. Therefore we come to the conclusion that, due to the singular behaviour which 
occurs at T = 0, we can not rely anymore on an expansion in powers of A when the temperature is lowered. It is 
conceivable that such an expansion could still be proper by accident for a specific phase, but it is unsafe for a general 
exploration of the various phases in competition. A possible partial cure for this problem could be to sum up the 
most divergent contributions, which are precisely those occuring in the FF phase. We have tried such an approach, 
but although it provides some improvement it clearly does not lead to a satisfactory situation. 

Therefore in an attempt to extend the simple approach around the TCP we will in the next section proceed to a 
Fourier expansion in the exact quasiclassical formulation of the problem. This will prove to be completely satisfactory. 



III. FOURIER EXPANSION 

We start from Eilcnberger's equations for the diagonal g(co, k, r) and off-diagonal f(ui, k, r) quasiclassical propaga- 
tors, which we simplify from the outset by taking h — 1 and m = 1/2. They read [11] : 

(lo + k.V)/(w, k, r) = A(r)g(w, k, r) 
(w-k.V)/ + (w,k,r) = A*(r) 5 (w,k,r) 
2k.V<?(o;,k,r) = A*(r)/(w,k,r) - A(r)/+(w,k,r) (7) 

where k is at the Fermi surface k — kp. Actually g is given in terms of / and / + by the normalization condition : 

g(u>, k, r) = (1 - f(u, k, r)/+ (w, k, r)) 1 / 2 (8) 

so the last equation results from the two first ones. These ones are also related [11] since : 

/*(-<*>, k,r) = / + (w,k,r) 5 *(-w,k,r) = -g(w,k,r) 

/*(w,-k,r) = / + (w,k,r) </*(w,-k,r)= ff (w,k,r) (9) 

In this paper we consider only an order parameter which varies only along the x axis. Accordingly / and g depend 
only on this variable. Moreover we assume that the order parameter is periodic in this direction which is the situation 
occuring in the FFLO transition. We restrict also ourselves to real order parameters since these have been found to 
correspond to the highest critical temperature in the vicinity of the TCP, and the LO solutions are also real, so this 
property is expected to be widely satisfied. Anyway the generalization to an intrinsically complex order parameter 
should not make much difficulties. 

Then we proceed to a Fourier expansion of this order parameter. Let us first assume, in order to present our method 
in the simplest case, that only the lowest harmonic is relevant. This amounts to take : 

A(x) = 2Acos(qx) (10) 

We will consider at the end of the paper the general situation, but we will actually find that, for our problem, the 
actual order parameter at the transition is very nearly a simple cosine. For fixed k Eilcnberger's equations are a set 
of first order differential equations for the variation of the Green's functions along k. So we take a reduced variable 
along this direction by setting r = kX, which gives k.V = d/dX and A(x) = 2Acos(QA) where we have introduced 
Q = kF<j cos 9 with 9 the angle between k and the x axis. Then we make a Fourier expansion of the Green's functions: 

f(X) = J2f n e inQX f + {X) = Y,ntte inQX g(X)^Y,9ne mQX (11) 

n n 

Explicit substitution of Eq.(ll) in Eilcnberger's equations Eq.(7) gives: 

A 

/ " = ^»Q (s "- 1+fl " +l) 

(fn-l + fn+l -/+_!-/+ +1 ) (12) 
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The solution of these equations have simple symmetry properties, which can be checked directly for example by 
generating explicitely the solution by a perturbation expansion. Actually they arise quite generally from the fact that 
we deal with an order parameter which is real and even (i.e. parity is not broken). This is more conveniently seen 
by taking the case where u> is real. However one has to keep in mind that we will deal finally with a complex u. 
Nevertheless the symmetry properties are still valid generally in this case. 

For real order parameter f(X), f + (X) and g{X) are real, which is consistent with Eilcnbcrger's equations. This 
implies /_„ = /*, f^ n = /+* and <?_„ = g* . Moreover for an even order parameter, Eq.(7) are unchanged when 

(k, r) is changed into (— k, — r) which shows that /, / + and g are also unchanged. Hence from Eq.(9) f(—X) = f + (X) 
and g(— X) = g(X), which leads finally to /+ = /_„ and g n — g_„. 

It is then convenient to make explicit the relation between /„ and /+ by introducing d n = (/„ — /+)/2i, which 
gives f n = (i — Lo/nQ)d n . We have then for the two quantities g n and d n (which are real for real lo) the recursion 
relations : 

nQA 

dn = 2~i 27^2^-1 + gn+l) 

A 

g n = —^(dn-i + d n+1 ) (13) 

It is clear from these equations that g n ^ only for even n, and d n ^ only for odd n, as it can be seen for example 
by generating the solution perturbatively. Moreover they satisfy g- n — g n and d- n = — d n . These equations are linear 
and must be supplemented by the normalization condition Eq.(8). The n = component is enough and it provides 
us precisely with the spatial integral go = J dr g s (ui, k, r) which we need in Eq.(l) to calculate the free energy : 

oo 

9t = 1 - E (9n9-n + fnflj = 1 " £ i^l + fn + f-n) (14) 
n=£0 n=l 

Now the interesting point is the large n behaviour of g n and d n . If for example we eliminate d n in Eq.(13) we obtain 
a linear recursion relation which links g n +2 to g n and g n -2- Since g- n = g n we have only to consider n > 0, but this 
becomes n > 2 when one takes into account that in Eq.(13) the relation for g is identically satisfied because d-\ = —d\. 
The general solution of such a recursion relation is a linear combination of two independent solutions. The large n 
behaviour is found from the recursion relation, which for A, \lu\ <g. \nQ\ simplifies into A 2 (g n+ 2 + gn-2) +n 2 Q 2 g n = 0. 
One sees that this equation has very rapidly growing solutions satisfying g n+ 2 S> g n 3> g n -2 an d behaving as 
<?2p+2 ~ (^1) P (2Q/A) 2p (p!) 2 . Naturally these solutions are not physically acceptable. On the other hand the recursion 
relation has also a solution satisfying g n +2 <C g n *C g n -2 and behaving as g2 P ~ (— l) p (A/2Q) 2p (l/p!) 2 , which is the 
physical solution we are looking for. This solution is found only if g and gi are related by a specific boundary 
condition. 

The very fast decrease of g n and d n provides an easy way to obtain a set of approximate solutions, which moreover 
converges rapidly to the exact one, all the more since these are <? 2 and d\ which come in Eq.(14) for the calculation 
of go- Since g n and d n are very small for large n we just take them to be zero beyond some fixed value. This serves 
as boundary condition. Then we work backward to obtain the whole set of Fourier components and normalize them 
properly through the normalization condition Eq.(14). Specifically we proceed as follows. Since the recursion relations 
are linear we rescale g n and d n in order to have convenient initial values. We set g n = goG n (which implies Go = 1) 
and d n = nQg D n and take as initial values G2N+2 — and D 2 n+i 7^ to be determined later. Then starting with 
p = TV we use for decreasing values of p the following recursion obtained from Eq.(13) : 

c 2 + (2p+l) 2 Q 2 

<^2p — — <-r2p+2 ^ ^p+l 

^-^i^^-i^^] as) 

down to Go- All the G's and D's are proportional to D 2 n+i, which is now found by enforcing Go = 1. Finally Eq.(14) 
gives explicitely for g : 

JV N 

g^ 2 = 1 + 2 G% + 2 J> 2 - (2P + l?Q 2 ]D% +1 (16) 
p— 1 p— 

When we let N — > 00 this equation provides the exact result for go- It is interesting to note that for these large n 
we have found that g n is proportional to A™. This makes a precise link between the expansion in powers of A we 
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discussed at the beginning and the Fourier expansion we are considering now. One can see our result as corresponding 
to resummations of infinite series, eliminating in this way the troubles mentionncd in section II occuring because 
coefficients in the Landau-Ginzburg expansion change sign as the temperature is lowered. One finds also that in the 
limit of large \u>\ ^> A, \nQ\, where one must recover the normal state Green's functions, one has g2 P ~ (— l) p (A/ u)) 2p . 
Naturally the recursion relations Eq.(15) are very convenient and very fast for a numerical implementation and in 
practice the situation is not very different from having an analytical expression for go . The only practical problem is 
linked to the determination of the square root in obtaining go from Eq.(16), but this is solved by noticing that, from 
the general spectral representation, one has Rego > when cj n > 0. 

The simplest of these approximations corresponds to take N = and it is given explicitely by : 

* = I' + ^^|^ (17, 

This is already a quite non trivial approximation. Since it is correct up to order A 2 it gives the proper location for 
the standard FFLO second order transition line. Moreover as we will see it gives qualitatively and semiquantitatively 
the correct results, with a first order transition down from the TCP which becomes a second order transition at low 
temperature in agreement with Ref. [10]. 

Although it is quite simple the calculation of the free energy has to be carried out numerically and naturally it is the 
same for all the higher order approximations. In practice it is convenient to make use in Eq.(l) of ln[T/T sp (p/T)] = 
1/N V — 7rT^sgn(w n )/a)„ to rewrite it as : 



N 



dr\A(r)\ 2 + AnT^J dcu Rc[< g (u, - tftk) > fc -1 + ^- _ ^_ )2 J rfr|A(r)| 2 ] (18) 



where we have made no assumption on the spatial dependence of A(r). In the present case Eq.(10) gives / (ir|A(r)| 2 = 
2A 2 . The form Eq.(18) is convenient for the frequency integration since the integrant behaves as u>~ for large ui, 
with a corresponding behaviour u>~ 3 in the Matsubara frequency summation. One may replace hi[T/T sp (p/T)] by 
ln[p/p sp (p/T)] where p sp (p/T) is the spinodal field for a given p/T, since p sp {p/T) /T sp (p/T) = p/T. Finally the 

angular average amounts to an integration over Q since < <7o(k) >&= / (d£lk /^)9o (k) = du go(Q = Jiqu) with 
u = cos9 and q = qkr/fi = Tiqkp I '(2mp). 

Since we are only interested in the transition from the normal to the FFLO state, in order to obtain the corresponding 
critical temperature T c we look for the highest temperature T (or effective field p.) at which f2 s — f2„ = 0. Precisely 
this leads to the following equation for T = T c : 



\n[T/T sp {fi/T)] = -Min— ^ / duj Rc[< g (uj - ifi, k) > fc -1 + , _ ._. 2 ] 



(19) 



Since for homogeneity the right-hand side of this equation is only a function of T/p,, A/p and q/p one has to minimize 
with respect to A/p and q/p, at fixed T/p. 

At T = this simplifies somewhat. The summation over Matsubara frequencies goes to an integral which is 
performed by a by parts integration. This gives for the critical field p : 

1 f°° - A 2 

ln[p/p sp (T = 0)] = -Min— J dujcjRc[< g (u) - ifi, k) > k -1 + ^_^ 2 ] (20) 

Finally we write down the self-consistency equation (or 'gap equation') for the order parameter, which comes out 
quite generally when the free energy is minimized with respect to variations of this order parameter. This equation 
gives the necessary feedback to Eq.(7) where A(r) can not be naturally an open function. In practice we will make 
little use of it since we will rather minimize directly the free energy with respect to the order parameter. The 
self-consistency equation [11] can be written as: 

oo . 

A n \n[T/T sp {p/T)} = 2nT ]T Rc[< f n (w - ip, k) > k -^] (21) 

by making the same transformation from T c0 to T sp (p/T) as we did to find Eq. (18). 
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IV. GENERAL ORDER PARAMETER 



We consider now the extension of our Fourier expansion, presented in section III, to a general order parameter: 

A(x) = A « e mQX ( 22 ) 

n 

We assume again a real order parameter, which implies A* = A_ n . We also restrict ourselves as before to an order 
parameter even with respect to x, which makes A„ real. Substitution as above in Eilcnberger's equations give the 
following generalization of Eq.(12): 

oc 

(u) + inQ)f n = (u>- inQ)f + = ^ A p (g n ^ p + g n+p ) 

OC 

2mQg n = ]T A p (/ n _ p + f n+p - f+_ p - f+ +p ) (23) 
P =i 

Here we have also assumed for the moment that the spatial average of the order parameter is zero, that is Ao = 0. 
Introducing again d n = (/„ — /+)/2i, we obtain the recursion relations: 

dn = ^ LU 2 +n 2Q2 H A p(5"-P + 9n+p) 
p=l 
_^ oc 

9n = —TjYl A p(^"-p + d «+p) (2 4 ) 

^ p=l 

We consider now the practical situation met in numerical use of these equations. In this case the number of Fourier 
components for the order parameter will be finite, so we have a maximum integer P for which A p = when p > P. 
We look again at the behaviour of the physical solution for d n and g n when n goes to infinity, and show that it is 
consistent with a fast factorial type decrease, as we have found in section III. Indeed in this case the dominant term 
in the sums found in the right hand side of Eq.(24) is the one where d n or g n has the smallest index n, which gives 
d n ~ —Apg n -p/(nQ) and g n <~ Apd n -p/(nQ). This leads to g n ~ (Ap/Q)™/ p /(n!) 1 / p . Although this is still a fast 
factorial type decrease, it gets slower when P increases. On the other hand the large n behaviour contains the power 
law dependence (Ap)"/ p , where in generic situations Ap is expected to be very small for large P. This is indeed what 
is found when one writes the self-consistency equation [11] Eq.(21) which gives A„ in terms of /„. This fast decrease 
of A„ corresponds to the standard situation, where there is no smaller physical length scale for A{x) than 1/q itself. 
However one can think of other particular situations, where this fast decrease of A„ does not hold and which should 
be dealt with specifically. Ultimately this convergence problem has to be handled numerically by making calculations 
for increasing P and looking when reasonable convergence has been achieved. This is what we will do below with our 
present problem of finding the location of the transition. 

Finally we make the same practical use of this fast convergence property as in section III. We take as boundary 
condition that g n and d n are zero beyond some fixed value N. This allows to calculate all the g n and d n within 
a common multiplicative factor, which is then found by the normalization condition Eq.(14). N is progressively 
increased until convergence has been obtained. The situation for solving the practical problem of finding the g„'s 
and the d„'s is less convenient than in section III. However we still have a linear problem for which very efficient 
numerical procedures are known. We have basically to handle a matrix. Instead of having a tridiagonal matrix, with 
just matrix elements right below and above the main diagonal we have now a band diagonal matrix with, in addition 
to the main diagonal, 2P diagonals with non zero matrix elements. 

To be more specific we have now to take into account that, in our problem, only odd Fourier components A2 P +i 
of the order parameter are nonzero. First we consider only order parameters with a zero spatial average Ao = 0, 
since taking a nonzero value amounts to mix in the order parameter of the uniform BCS phase, which is energetically 
unfavorable. Hence it is reasonable to assume that similarly a nonzero A is unfavorable. Next we can see for example 
by a iterative treatment to all orders, in order to obtain an exact solution of Eilenberger's equations, that we have 
only odd components. Indeed if we start with the simple A(x) = 2Acos(qx) that we considered in section III, we 
generate only odd Fourier components in f(X) and even components in g(X) as we have seen. Now this f(X) in turn 
generates only odd components for A(x) from the self-consistency equation Eq.(21). But from Eq.(24) this is again 
completely compatible with only odd components for f(X) and even for g(X). Naturally it can also be seen directly 



7 



from the starting Eilcnberger's equations that such a solution is a consistent one. We note that such a solution with 
odd components means that, by shifting x by 7r/(2g), we obtain an order parameter which is odd with respect to x, 
in the same way as it transforms Eq.(lO) into 2Asin(qx). 

Then it results from Eq.(24) that, just as in section III, g n ^ only for even n, and d n ^ only for odd n. In the 
same way we set g n = goG n (implying Go = 1) and d n = nQg D n , and we have again G_„ = G„ and = D n . It 
is now convenient to include g n and d n into a single unknown vector V n , defined by V 2p = G 2p and V2 P +\ = £> 2p +i. 
Then Eq.(24) can be merely written as M mn V n = A m with A n = — A„ and the matrix M given by: 

M 2n .2n = 2n 
M 2n+ i,2n+i =uj 2 + (2n + 1) 2 Q 2 

M2 m +l,2n = A 2 ( m+ „) + 1 + ^|2(m-n)+l| 

M 2m ,2n+i = (2n + l)[A 2 ( m+n )+i - A| 2 („_ m )+i|] (25) 

with m,n > 1. 

As explained above we truncate the infinite matrix M by to, n < A max , which gives a matrix with finite order iV max . 
The corresponding linear equation for V n , with n < A max , can be solved numerically by efficient standard routines, 
since as mentionned above the matrix M has a generalized band diagonal form. Once this is done, go is still obtained 
from Eq.(16), the free energy calculated from Eq.(18) and the critical temperature obtained by minimization. Finally 
the whole procedure is repeated for increasing values of A max until convergence is achieved. In the next section we 
will display the corresponding numerical results. 

V. NUMERICAL RESULTS 

We present now the results of our numerical implementation of the above procedure. In the first susbsection below 
we restrict ourselves to an order parameter with only the lowest harmonic as it is given by Eq.(lO). The general case 
is considered afterwards. 



A. Lowest harmonic 

We first give in Fig. 1 the results for the calculation of the critical temperature, down from the TCP. Rather than 
giving T c (p), we cover for convenience the (/2, T) plane in polar coordinates, rather than cartesian coordinates, and 
give the critical temperature T c (p,/T) as a function of ji/T, equivalent to a polar angle. More precisely we plot its 
ratio T c /Tfflo to the FFLO critical temperature Tfflo{p-/T) obtained for the same value of the ratio ji/T. This 
is more suited to the present situation since we find this ratio to be always near unity. However to make the graph 
more readable, we give on the x-axis, instead of ft/T, the value of TpFLoifi/T) itself, compared to the standard BCS 
critical temperature T c0 , found for fi, = (this corresponds to go along the standard FFLO transition line). Naturally 
when our result for T c /Tfflo goes below 1, this means that the first order transition is less favorable than the second 
order one, so when the temperature is lowered there is actually a switch from first to second order when one finds 
that T c /T FFLO = 1. 

We give the results of the calculation with increasing values of N max going up to 6. The approximation iV max = 1 
corresponds to the explicit result Eq.(17) for g . As already mentionned above, it is correct up to second order in 
A and consequently it gives the correct location for the FFLO transition. Moreover we see that it gives already the 
proper result semi-quantitatively for the order of the transition, since it is gives a switch from first to second order 
when the temperature goes below Tfflo/T c q = 0.195. The next approximation A max = 3 for odd A max is already 
quite good quantitatively since it gives 0.063 for the above ratio. Full convergence is obtained for A max = 5 where we 
find Tfflo/Tco = 0.076 in very good agreement with Matsuo et al [10]. For completeness we give also in Fig. 1 our 
results for even A max . It is less natural, from the structure of the recursion equations, to truncate them in this way. 
Hence it is not so surprising that the approximation A^ max = 2 is much worse than iV max = 1 since it does not even give 
a switch from first to second order for the transition. Nevertheless we have naturally convergence when we increase 
A max , and indeed we find that A max = 4 is already very good since the switch is located at Tfflo/T c o = 0.074, while 
Am ax = 6 is completely converged. 
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FIG. 1. Critical temperature T c (p,/T) for the transition to the order parameter given by Eq.(10), divided by the FFLO critical 
temperature Tfflo(P-/T) obtained for the same value of ft/T. On the x axis, instead of fi/T, we have given Tfflo{U) /T c q, 
where T c o is the maximal critical temperature, obtained for fl — 0. 

A noticeable feature of Fig. 1 is that the ratio T c /Tpp^o stays always very near unity, while one would have 
expected to find it larger since there is no obvious relation between the order parameters of the first and the second 
order transitions. This behaviour is also found [7] near the TCP. A natural conclusion from this feature is to say 
that the first order transition is actually always quite near to be a second order one. We can check in our results if 
this interpretation is a coherent one by looking at the size of the order parameter (more precisely its maximal value 
as a function of spatial position), that is essentially the value of Ai at the first order transition. If the first order 
transition is nearly a second order one, it should be small compared to a typical gap Ao in the uniform BCS phase. 
Equivalently, in Fig. 2, we compare it to \x since it is of order Ao in all the range we are interested in (at T — the 
FFLO result is /2 = 0.754Ao). Our results show clearly that the size of the order parameter at the transition is of 
order of the one deep in the standard BCS phase, so it is not at all possible to consider that the first order transition 
is nearly a second order one. 

0.6 
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0.4 



0.2 
0.1 


0.1 0.2 0.3 0.4 0.5 0.6 

rFFLo/r c o 

FIG. 2. Full line: size of the component Ai of the order parameter for the first order transition, compared to p,, when all 
higher order harmonics are taken equal to zero. Dashed line: same quantity for an order parameter where A3 is also non zero. 

Finally it is also of interest to compare the reduced wavevectors of the order parameter for our first order transition 
and for the standard second order FFLO transition. This is done in Fig. 3 where it is seen that, although there are 
differences, they not very large so that the reduced wavevector is rather similar for the two transitions, in contrast 
with the size of the order parameter. 
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FIG. 3. Full line: optimal reduced wavevector q of the order parameter at the first order transition, for the converged solution, 
as a function of temperature, when only the component Ai is different from zero. Long dashed line: same result when A3 is 
also non zero. Short dashed line: corresponding result for the second order FFLO transition 



Naturally it is not consistent to keep only the lowest harmonic in the order parameter, as it is immediately seen 
from the self-consistency equation (21). Hence we consider now the effect of higher harmonics. In a first step we have 
explored the effect at the transition of the inclusion of A 3 . We have found it quite small. This would imply normally 
to stop the exploration at this stage since one expects the effect of harmonics A5 and higher to be even smaller. 
However one might wonder whether this result is not somewhat accidental [20] and specific to A3. This is especially 
a concern in the vicinity of the switching temperature from first to second order, where A3 is particularly small (see 
Fig. 5). It could be that higher order harmonics dominate in this region leading to a quantitative change of the first 
order transition line. Hence, in order to eliminate any doubt on this question, we have explored the effect of taking 
A5 and A7 to be different from zero. Our study shows that these harmonics give also a very small contribution. 
Hence our conclusion is that the optimal order parameter remains always very close to a simple cosine in the whole 
temperature range from the tricritical point down to zero temperature. 

Our numerical procedure is to use directly the free energy Eq.(18) by taking as an ansatz our form for the order 
parameter, with either three or four odd Fourier components. More precisely we maximize, with respect to q, Ai with 
i = 1, 3, 5, 7, the critical temperature from the generalization to our case of Eq.(19), as explained in section III. We 
then check that our optimal form satisfies the gap equation. We have also performed calculations by making use only 
of the gap equation. The results are not significantly different from the ones we display below, and most of the time 
agree with them within numerical accuracy . From a practical point of view, we have chosen high enough values of 
N max , typically N max = 12, so that numerical results are insensitive to changes in iV max . 

We give first in Fig. 4 our result for the critical temperature of the first order transition. The effect of all our higher 
order harmonics can be only barely seen in the figure, as compared to our calculation with only the lowest harmonics 
Ai, already given in Fig. 1. 



B. Higher harmonics 
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FIG. 4. Dashed line: critical temperature for the first order transition for a one dimensional order parameter form with four 
odd Fourier component Ai with i = 1, 3, 5 and 7. Full line: same result for the simple cosine ansatz, where only Ai is different 
from zero. 

Next we display in Fig. 5, as a function of temperature, the values of the higher order harmonics A 3 ,A 5 and A7 for 
the optimal order parameter. It is seen that they are always quite small compared to Ai. Nevertheless, around and 
below the switching temperature, A3 and A5 are of the same order while one would have rather expected A5 to be 
small compared to A 3 (note that anyway these results are physically irrelevant below the switching temperature since 
they are for the first order transition, while the actual transition is second order). On the other hand A7 is always 
negligible compared to A3 and A5, except near the TCP where anyway A5 and A7 are essentially zero. Finally Fig. 
5 shows also the results for the optimal A3 and A5 when A7 = 0. The difference with the preceding results is within 
numerical error. Similarly our result for A3 when A5 = A7 = (not shown) are also essentially indistinguishable 
from the result displayed in Fig. 5. 
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FIG. 5. Optimal values for the amplitudes A3, A5 and A7 in the order parameter as a function of temperature (full lines). 
For comparison the results for A3 and A5 when A7 = are also given as dashed lines. 

VI. CONCLUSION 

In this paper we have shown that performing a Fourier expansion in the quasiclassical Eilenberger's equations 
provides a very efficient way to study the transition from the normal state to the FFLO phases in 3 dimensions, 
at least in the vicinity of the transition. We have applied this technique to the case of the transition to the one- 
dimensional 'planar' order parameter and we have found perfect agreement with the earlier work of Matsuo et al. In 
particular we have rederived their remarkable result that, when the temperature is lowered, the transition switches 
from first to second order. We have shown in detail that, in the case of the first order transition, the order parameter 
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is nevertheless dominated by its lowest order Fourier component, in somewhat surprising contrast to what one might 
guess for such a transition. This feature contributes naturally to make our Fourier expansion very rapidly converging. 

However the strength of our method is not so much displayed in this case of a one-dimensional order parameter. 
Its major interest is rather that our approach can be fairly easily generalized to more complex order parameters, 
with full three-dimensional spatial dependence. As shown by Larkin and Ovchinnikov these order parameters come 
in competition and, in the case of a first order transition, it is not clear that they are not more advantageous than the 
standard second order FFLO phase transition. We will indeed show, in forthcoming work, that this is the case at low 
temperature in 3 dimensions. Finally another interest of our approach is to provide some insight, even if approximate, 
in the analytical structure of the theory, as we have seen by providing explicit approximate analytical solutions. We 
believe that this might be helpful in a theoretical situation where the intrinsic non linearity of the equations forces 
mostly to purely numerical work. 

* Laboratoire associe au Centre National de la Recherche Scientifique et aux Universites Paris 6 et Paris 7. 
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